In [1]:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 5/4/2025
@author: ifenty
"""
import numpy as np
from pathlib import Path
import xarray as xr
import matplotlib.pyplot as plt;
In [3]:
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
In [5]:
def calc_adxx_length(adxx):
tmp = []
# number of time levels
for i in range(adxx.shape[0]):
tmpv = np.inner(np.ravel(adxx[i]), np.ravel(adxx[i]))
tmp.append(tmpv)
return np.array(tmp)
In [6]:
def plot_adxx_length(adxx, lags, t=''):
plt.figure()
adxx_len = calc_adxx_length(adxx)
plt.plot(lags,adxx_len,'k.-')
plt.title('|adxx ' + t + '| vs. lag')
In [7]:
def plot_gradient_field(adxx_DA, central_longitude=-180):
field_to_plot_max = np.nanmax(np.abs(np.ravel(adxx_DA)))*0.6
# Plot
fig = plt.figure(figsize=(15, 5))
ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=central_longitude))
pcm = ax.pcolormesh(adxx_DA.lon, adxx_DA.lat, adxx_DA,
vmin = -field_to_plot_max,
vmax = field_to_plot_max,
transform=ccrs.PlateCarree(), cmap="RdBu_r")
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.colorbar(pcm, orientation='vertical')
plt.title(adxx_DA.name + ' at lag ' + str(adxx_DA.lag.values))
plt.show()
Define Paths¶
Paths of EMU gradient experiments
In [8]:
emu_output_basedir= Path('/home/ifenty/incoming/gradients')
gradient_dirs = ['emu_adj_48_48_1_15_743_1','emu_adj_48_48_1_2_209_1','emu_adj_48_48_1_56_1069_1','emu_adj_48_48_1_56_1072_1']
Find all EMU adjoint experiment directories¶
In [11]:
gd_output_dirs = []
gd_files = dict()
for gd in gradient_dirs:
tmp = emu_output_basedir / gd / 'output'
print(tmp)
gd_output_dirs.append(tmp)
gd_files[gd] = np.sort(list(Path(tmp).glob('*nc')))
print(gd_files[gd])
adxx_files = ['adxx_qnet','adxx_tauu','adxx_tauv']
/home/ifenty/incoming/gradients/emu_adj_48_48_1_15_743_1/output
[PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_15_743_1/output/adxx_qnet.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_15_743_1/output/adxx_tauu.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_15_743_1/output/adxx_tauv.nc')]
/home/ifenty/incoming/gradients/emu_adj_48_48_1_2_209_1/output
[PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_2_209_1/output/adxx_qnet.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_2_209_1/output/adxx_tauu.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_2_209_1/output/adxx_tauv.nc')]
/home/ifenty/incoming/gradients/emu_adj_48_48_1_56_1069_1/output
[PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_56_1069_1/output/adxx_qnet.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_56_1069_1/output/adxx_tauu.nc')
PosixPath('/home/ifenty/incoming/gradients/emu_adj_48_48_1_56_1069_1/output/adxx_tauv.nc')]
/home/ifenty/incoming/gradients/emu_adj_48_48_1_56_1072_1/output
[]
In [31]:
# loop through each EMU adjoint gradient experiment directory
# directory 3 and 4 are almost identical, so skip 4
# choose an EMU experiment to process
for dd in range(3):
# for plotting, select the central longitude
if dd == 0: # tropical pacific experiment
central_longitude = -180 # degrees east
else: # Nantucket and N. Atlantic experiments
central_longitude = -45 # degrees east
# make maps of the gradients for several lag
lags_to_show = [-1, -6, -36, -72]
adxx_qnet_DA = xr.open_dataset(gd_output_dirs[dd] / 'adxx_qnet.nc')['adxx_qnet']
adxx_tauu_DA = xr.open_dataset(gd_output_dirs[dd] / 'adxx_tauu.nc')['adxx_tauu']
adxx_tauv_DA = xr.open_dataset(gd_output_dirs[dd] / 'adxx_tauv.nc')['adxx_tauv']
# plot gradients for qnet, tauu, tauv
for adxx_DA in [adxx_tauu_DA, adxx_tauv_DA, adxx_qnet_DA]:
# plot different lags
for ll in lags_to_show:
plot_gradient_field(adxx_DA[ll], central_longitude=central_longitude)
# plot the length of the gradient vector (crudely calculated)
tmp = calc_adxx_length(np.where(np.isnan(adxx_DA.values),0, adxx_DA.values))
plot_adxx_length(tmp[-72:], adxx_DA.lag[-72:], adxx_DA.name);plt.grid()
plt.xlabel('days')
plt.ylabel('gradient length')
# close the DataArray file
adxx_qnet_DA.close()
adxx_tauu_DA.close()
adxx_tauv_DA.close()
In [202]: